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ABSTRACT 

The  main  purpose  of  this  work  is  to  give  explicit  sparsity¬ 
preserving  SOR  (successive  overrelaxation)  algorithms  for  the  solution 
of  separable  quadratic  and  linear  programming  problems.  The  principal 
and  computationally  distinguishing  feature  of  the  present  SOR  algorithms 
is  that  they  preserve  the  sparsity  structure  of  the  problem  and  do  not 
require  the  computation  of  the  product  of  the  constraint  matrix  by  its 
transpose  as  is  the  case  in  earlier  SOR  algorithms  for  linear  and 
quadratic  programming. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  important  optimization  problems  have  linear  constraints  and 
objective  functions  which  are  quadratic  or  linear.  Very  often  these 
problems  are  very  large  but  sparse.  Conventional  methods  such  as  the 
simplex  method  and  other  pivotal  methods  may  not  be  able  to  handle 
such  problems  because  of  their  size  and  because  the  sparsity  structure 
of  the  problem  may  be  quickly  lost  when  these  methods  are  used.  We 
propose  here  a  different  class  of  methods,  successive  overrelaxation 
(SOR)  methods,  which  can  handle  large  problems  while  preserving  their 
sparsity.  SOR  methods  have  been  widely  and  successfully  used  in  the 
solution  of  linear  systems  of  equations,  but  rarely  in  the  solution  of 
optimization  problems. 


The  responsibility  for  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


SPARSITY-PRESERVING  SOR  ALGORITHMS 
FOR  SEPARABLE  QUADRATIC  AND  LINEAR  PROGRAMMING 

0.  L.  Mangasarian 

1.  Introduction 

Recently  iterative  SOR  methods  have  received  widespread  attention 
in  the  solution  of  the  symmetric  and  nonsymmetric  linear  complementarity 
problem  [4,6,3,11,15,16,1],  quadratic  and  linear  programming  problems 
[12,13].  In  the  case  of  the  latter  two  problems  which  are  our  principal 
concerns  here,  the  recently  proposed  SOR  algorithms  do  not  preserve  any 
sparsity  that  the  original  problems  may  have  had.  This  is  due  to  the 
fact  that  algorithms  as  presented  in  [12]  require  the  product  of  the 
constraint  matrix  by  its  transpose,  which  can  cause  loss  of  both  sparsity 
and  accuracy.  In  this  work  we  shall  present  some  explicit  realizations 
of  the  algorithms  of  [12,13]  which  will  not  require  the  multiplication 
of  the  constraint  matrix  by  its  transpose.  These  computationally 
improved  realizations  which  follow  from  the  algorithms  of  [12]  have  not 
been  given  explicitly  before.  The  absence  of  such  sparsity-preserving 
algorithms  has  been  a  critical  factor  in  preventing  the  application  of 
SOR  methods  to  many  large  important  but  highly  structured  problems  such 
as  economic  equilibrium  problems,  transportation  and  network  flow 
problems.  In  addition  some  of  the  present  realizations  of  the  SOR 
algorithms  (e.g.  (14)  and  (32)  below)  require  only  simple  operations  on 
the  rows  of  the  constraint  matrix,  and  hence  very  large  problems  can  be 
tackled  by  such  SOR  realizations,  because  only  linear  row  arrays  are 
needed  in  the  computations.  These  advantages  become  even  more  pronounced 
if  these  linear  row  arrays  are  sparse  and  hence  can  be  stored  in  packed 


form. 


The  paper  is  organized  as  follows.  In  Section  2  we  give  an  SOR 

algorithm  for  the  symmetric  linear  complementarity  problem  or  equivalently 
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for  the  quadratic  programming  problem  with  nonnegativity  constraints  only. 
This  is  a  special  case  of  the  general  algorithm  presented  in  [11]  but 
given  here,  in  a  simple  explicit  form  in  terms  of  the  rows  of  the  matrix 
defining  the  problem,  principally  to  make  it  preserve  problem  sparsity. 

In  Section  3  we  consider  a  separable  quadratic  programming  problem  and 
give  a  version  of  the  SOR  algorithm  of  [12]  which  does  not  require  multi¬ 
plication  of  the  constraint  matrix  by  its  transpose.  Hence  this  present 
form  of  the  algorithm  is  now  ideally  suited  for  large  sparse  problems.  In 
Section  4  two  sparsity-preserving  SOR  algorithms  for  linear  programming 
are  given.  One  is  based  on  finding  the  ''smallest"  optimal  primal-dual 
solution  (LPS0R1)  [13]  and  the  other  is  based  on  perturbing  a  linear 
program  to  a  separable  quadratic  program  and  then  solve  the  latter  by  the 
method  of  Section  3  (LPS0R2)  [12],  Computational  experience  with  a 
version  of  LPS0R1  [9]  on  problems  of  size  up  to  800  constraints  and  1000 
variables  and  the  non-sparsity-preserving  version  of  LPS0R2  [12]  have 
been  very  encouraging.  It  is  hoped  that  further  refinements  will  make 
SOR  methods  simple,  robust  and  commercially  viable  methods  for  solving 
very  large  separable  quadratic  and  linear  programs. 

We  briefly  describe  now  the  notation  used.  All  matrices  and  vectors 
are  real.  For  the  mxn  matrix  A  we  write  AeRm*n  and  denote  row  i 
by  Aj,  column  j  by  A*j  and  the  element  in  row  i  and  column  j  by 
A^j.  For  x  in  the  real  n-dimensional  Euclidean  space  Rn,  element  i 
is  denoted  by  x^  and  x+  will  denote  the  vector  with  components 
(x+).  =  max  {x^,0},  ial,...,n.  All  vectors  are  column  vectors  unless 
transposed  by  the  superscript  T.  ||x||  will  denote  the  2-norm, 
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(x^x)*5  =  (  y  x?)*5.  A  matrix  C  In  Rnxn  is  positive  semidefinite  if 
j=l  J 

x^Cx  >  0  for  all  x  in  Rn  and  positive  definite  if  x^Cx  >  0  for 
all  nonzero  x  in  Rn.  For  brevity  we  shall  sometimes  omit  mentioning 
the  dimensionality  of  a  vector  or  matrix,  it  being  obvious  from  the 
context.  The  vector  e  will  be  a  vector  ones  in  a  Euclidean  space  of 
appropriate  dimension.  For  a  twice  differentiable  function  4>:RmxRn-*R, 
Vu4>(u,v)  will  denote  the  mxl  gradient  vector  with  elements 

3<%-v )■,  i=l,...,m,  Vj(u,v)  will  denote  the  nxl  gradient  vector 

oil  *  V 


with  elements  ,  i=l,...,n,  V<j>(u,v)  = 


Vu<Ku,vJ] 

VY$(u,v) 


and 


V^u.v) 


will  denote  the  Hessian  in  R(n+m)x(n+m)  wi th  submatrix  components 
denoted  as  follows 


V2<|>(u,v) 


Vuu4>(u,v)  V^u.vf 
\u<j»(u,v)  Vyv<j)(u,v) 
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2.  SOR  Algorithm  for  the  Symmetric  linear  Complementarity  Problem 
We  consider  here  the  problem  of  finding  z  In  Rn  such  that 

Mz  +  q  >  0,  z  >  0,  zT(Mz+q)  =0  (1) 

where  M  is  a  symmetric  matrix  in  R^xk  and  qeR^.  Conditions  (1) 
are  [10]  the  necessary  optimality  conditions  for  the  quadratic  program¬ 
ming  problem 

minimize  izTMz  +  qTz  subject  to  z  >  0  (2) 

z*Rk  2 

Conditions  (1)  are  sufficient  for  z  to  solve  (2)  whenever  M  is 
positive  semidefinite  [10]. 

In  [4,6,11]  iterative  SOR  methods  have  been  proposed  for  solving 
(1),  but  without  paying  any  special  attention  to  possible  sparsity  that 
the  problem  may  have.  We  give  below  an  SOR  algorithm  based  on  that  of 
til]  in  which  any  sparsity  that  exists  is  left  undisturbed.  If  we 
define 


0(z):=  |-zTMz  +  qTz 


(3) 


then  the  SOR  algorithm  for  solving  (2)  can  be  represented  as  a  gradient 
projection  algorithm  of  the  following  type 


*  (zj-(l){72e(zi))Tj  vz  6(z]+1, 

<3 


zi+1  z1  z1!! 

‘  j-i,zj . *V'+ 

j— !»•»»«, k 


(4) 


where  u  is  the  relaxation  factor  or  stepsize  that  must  be  in  the  open 
interval  (0,2)  and  i  represents  the  ith  iteration.  More  specifically 
we  have  the  following. 
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LCPSOR  Algorithm 


Choose  zu«Rj,  w£{0,2).  Having  z1  compute  z^  as  follows: 


z1+1  =  (z1-wMt!(  M.0zj+1+  jM.#z!+q.)) 

J  '  J  ^  JA  £  JA  A  V'+ 
for  J>1  w 

J"  I  S  •  «  •  «  |N 


(5) 


The  following  convergence  theorem  follows  directly  from  [11]. 

Theorem  1:  LCPSOR  Convergence 

(i)  Let  M  be  symmetric.  Each  accumulation  point  of  (5)  solves 
(1).  If  in  addition  M  is  positive  semidefinite  then  each 
accumulation  point  of  (5)  solves  (2)  as  well. 

(ii)  Let  M  be  syrmetric  and  positive  semidefinite  and  such  that 

Mz  +  q  >  0  for  some  z  e  Rn  (6) 

Then  the  sequence  {z*}  of  the  LCPSOR  algorithm  (5)  is 
bounded  and  has  an  accumulation  point  that  solves  both  (1) 
and  (2). 

(iii)  Let  M  be  symmetric  and  positive  semidefinite  and  such  that 
problem  (1)  (or  equivalently  problem  (2))  has  a  nonempty 
bounded  solution  set.  Then  the  sequence  {z1}  of  the  LPSOR 
algorithm  (5)  is  bounded  and  has  an  accumulation  point  that 
solves  both  (1)  and  (2). 

(iv)  Let  M  be  symmetric  and  positive  definite  then  the  sequence 
(z1)  of  the  LCPSOR  algorithm  (5)  converges  to  the  unique 
solution  5  of  (1)  and  (2). 
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Proof 

Parts  (1),  (11)  and  (1 v)  follow  from  Theorem  2.1,  Theorem  2.2  and 
Corollary  2.2  of  [11]  respectively.  To  establish  (ill)  we  note  that 
from  Lenina  2.3(b)  of  [11]  that  if  the  sequence  {z^}  of  (5)  is 

_  I. 

unbouned  then  there  exists  a  yeR  such  that 

0  i  y  >  O,  My  *  0,  qTy  <  0 

This  contradicts  the  boundedness  assumption  on  the  solution  set  of  (1) 
since  if  z  solves  (1)  then  z  +  Xy  also  solves  (1)  for  all  X  >  0 
because  I  +  Xy  0,  M(z+Xy)  +  q  >  0  and 

0  <_  (z+Xy)T(M(z+Xy)+q)  =  Xq^y  ^  0.  □ 
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3.  SOR  Algorithm  for  Separable  Quadratic  Programming 
We  consider  here  the  separable  quadratic  program 

minimize  ixTDx  +  cTx  subject  to  Ax  >  b,  x  >  0  (7) 

xcRn  2 

where  D  is  a  positive  diagonal  matrix  in  Rnxn,  AeRmxn,  ccRn  and 
bcRm.  For  more  general  quadratic  programs  see  [12].  Associated  with 
this  quadratic  program  is  the  dual  quadratic  program  [5,17,10] 

maximize  -i-xTDx+bTu  subject  to  Dx  -  ATu- v  +  c  =  0,  (u,v)  >  0  (8) 

(x,u,v)eRn+m+n  2 

which  upon  elimination  of  x  by  using  the  constraint  relation 

x  =  D_1(ATu+v-c)  (9) 

gives 

minimize  i(ATu+v-c)T0”'*  (ATu+v-c)  -bTu  subject  to  (u,v)>0  (10) 

(u.vJcR^  2 

This  problem  (10)  is  now  precisely  of  the  form  (2)  and  the  LCPSOR 
algorithm  (5)  can  be  applied  to  it  easily.  Because  our  principal 
interest  here  is  sparsity  preserving  we  shall  spell  out  the  algorithm 
for  solving  (10)  explicitly.  Define  the  objective  function  of  (10)  as 

4>(u,v):=  y(A^u+v-c)^D"^  (A^u+v-c)  -  b*u  (11) 

then 

AD-^ (A^u+v-c)  -  b 
D_1(ATu+v-c) 


V$(u,v)  = 


(12) 
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Remark  2 

The  iteration  (14)  is  very  well  suited  for  matrices  A  which  have 
a  pronounced  row  structure,  for  example  if  A  is  sparse  and  the  nonzero 
elements  of  each  row  can  be  easily  located  without  search.  On  the 
other  hand  if  the  matrix  A  has  a  pronounced  column  structure,  then 
the  following  alternate  but  equivalent  iteration  to  (14)  may  be 
preferrable: 


,T 


. . >v-cvv,.rv,t 

J 

j”l>***# ,m 


(14f) 


'j+1  =  (Vj-w(u1+1  A.j+Vj-Cj))+  »  J=1 . n 


Theorem  2:  QPSOR  Convergence 

•  • 

(i)  Each  accumulation  point  (u,v)  of  the  sequence  {(u^v1)} 
generated  by  the  QPSOR  algorithm  (14)  solves  (10),  and  the 
corresponding  x  determined  by  (9)  solves  the  quadratic 
program  (7). 

(ii)  Let  the  feasible  region  of  the  quadratic  program  (7)  satisfy 
the  Slater  constraint  qualification 

{x]Ax>b,  x>0}  f  0  (15) 

1  • 

Then  the  sequence  { (u1 , v1 ) >  of  the  QPSOR  algorithm  (14)  is 
bounded  and  has  an  accumulation  point  (u,v)  and  the 
corresponding  x  determined  by  (9)  solves  (7). 


1 


Proof 


(1)  Follows  from  Theorem  1  (i )  and  the  duality  theory  of  quadratic 
programming  [10,  Theorem  8.2.5], 

(11)  Because  of  (15)  there  exist  a  6  >  0  such  that  the  perturbed 
positive  definite  quadratic  program 

minimize  ix^Dx  +  cTx  subject  to  Ax  >  b  +  ed,  x  >  ed 
xcRn  *  ~ 

has  a  solution  x*Rn  with  corresponding  multipliers 
(u,v)cRm+n  that  satisfy  the  Karush-Kuhn-Tucker  conditions 

DX  +  c  -  ATD  -  v  =  0,  Ax  >  b  +  ed,  x  >  ed,  £J  >  0,  ?  >  0 
u^(Ax-b-ed)  =  0,  vTx  =  0 


Hence 


X  =  0"1(ATG+v?-c)  >  ed  >  0 
AD“^(ATu+v-c)-b  >  ed  >  0 


Conditions  (16)  are  equivalent  to  condition  (6)  for  problem  (10).  Hence 
by  Theorem  l(ii)  the  sequence  {(u^.v1)}  of  the  QPSOR  algorithm  (14)  is 
bounded  and  has  an  accumulation  point  (u,v)  which  solves  (10).  Hence 
the  corresponding  x  determined  by  (9)  solves  (7).  □ 
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4.  SOR  Algorithm  for  Linear  Programming 

We  consider  finally  the  dual  linear  programs 

.  T 

minimize  c'x  subject  to  Ax  >  b,  x  >  0  (17) 

xcRn 

and 

maximize  b^u  subject  to  A^u  <  c,  u  >  0  (18) 

u«Rm 

where  AeRm><n,  ceRn  and  beRm.  It  is  well  known  [2]  that  solving 
either  (17)  or  (18)  is  equivalent  to  solving  both  (17)  and  (18)  which 
in  turn  is  equivalent  to  solving  the  linear  complementarity  problem 

Ny  +  p  >  0,  y  >  0,  yT(Ny+p)  »  yTp  =  0  (19) 

where 


k  =  n  +  m 


(20) 


Note  that  N  is  skew  symmetric,  that  is  N  +  NT  *  0  and  hence 
yTNy  =  0.  As  proposed  in  [13]  one  way  of  solving  the  linear  program 
(17)  is  to  find  the  closest  point  to  the  origin,  in  the  2-norm,  of  the 
solution  set  of  (19).  That  is  we  shall  solve  the  quadratic  program 


Minimize  -i  J(y li ^  subject  to  Ny  +  p  >_  0,  y  £  0,  p^y  <  0  (21) 

y«Rk 

Note  that  under  the  constraints  Ny  +  p  >  0,  y  >  0,  the  constraint 
p^y  <  0  is  equivalent  to  pV  a  0  since  0  >  ■  y^(Ny+p)  >  0. 

The  dual  to  the  quadratic  problem  (21)  is  [5,17,10] 
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maximize  -  2"}|y  H2-pTs  subject  to  y  -  NTs  -  t  +  Bp  ■  0,  (s,t,B)>0  (21) 

(y,s,t,8)eR3k+1 

Elimination  of  y  by  using  the  constraint  relation 

(j]  -  y  -  NTs  +  t  -  BP  (22) 

gives  the  quadratic  program 

minimize  i JjNTs-8p+tj}2  +  pTs  subject  to  (s,t,8)>0  (23) 

(s,t.6)£R2k+1 

Problem  (23)  and  consequently  problem  (21)  can  be  solved  by  the  SOR 
method  of  Section  3.  For  that  purpose  it  is  convenient  to  let 
<J>(s,t,B)  equal  the  objective  function  of  (23)  that  is 

$(s,t,B):=  ]rllNTs-BP+t!|2  +  PTs  (24) 

and  consequently 

V<J>(s,t,8) 


N(NTs-Bp+t)  +  p 

NTs  -  Bp  +  t  (25) 

-pT(NT$-BP+t) 
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It  is  obvious  from  (26)  that  7  <j>(s,t,B)  is  positive  semidefinite. 
We  can  now  state  an  SOR  algorithm  for  solving  (23)  based  on  QPSOR. 


LPS0R1  Algorithm 

r t _ _  /.0  4.0  o0\  n2k+l 


Choose  (s°,t°,B°)  cR^+  ,  we  (0,2).  Having  (s’.t’.B1)  compute 
(s1+W)  as  follows: 


»  (H.(  Jj'  (NT).  s|+1+  t  (NT)..s1-bW)+P,))+.  0=1 ... 

J  J  IlNjll2  J  1=1  4  ‘  44  J  + 


for  j>l 


t'+1  -  (t'-U(NTs1+,-8ip+ti))+  (27) 

61+1  .  (gl+  “  PT(»T*1+,-8Vt1+1» 
llPlI2 


Parts  (i)  and  (ii)  of  the  following  convergence  theorem  follow 
directly  from  Theorem  2(i)  and  Theorem  l(ii)  above  respectively. 


Theorem  3:  LPS0R1  Convergence 


(i)  Each  accumulation  point  (s,t,B)  of  the  sequence 

{(s^.t1^)}  generated  by  the  LPS0R1  algorithm  solves  the 

dual  program  (23)  and  the  corresponding  (*)  determined 

by  (22)  solve  the  dual  linear  programs  (17)-(18). 

(ii)  If  there  exist  (s,t,B)  e  R2k+1  satisfying  7<j»(s,t,B)  >  0, 

•  •  • 

then  the  sequence  { (s1 , t1 , B1 ) )  generated  by  the  LPS0R1 


algorithm  is  bounded  and  has  an  accumulation  point. 


I 


I 


•  i 


3 


^  -I 
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Note  that  Theorem  2(ii)  does  not  apply  here  because 
Ny  +  p  >  0,  y  >  0  imply  that  p^y  >  0  and  hence  we  cannot  satisfy 
the  Slater  constraint  qualification  that  there  must  exist  a  y 
satisfying  Ny  +  p>0,y>0  and  p^y  <0.  We  further  note  that 
the  condition  74>(s,t,g)  >  0  is  sufficient  but  not  necessary  for 
the  boundedness  of  the  sequence  {(s^,t*,S^)}.  Numerical  experiments 
have  revealed  no  serious  problems  with  unboundedness  of  the  sequence 
{(s1 .t1 .B1 )}  generated  by  the  LPS0R1  algorithm. 

We  conclude  by  giving  a  sparsity-preserving  version  of  the  SOR 
algorithm  for  solving  a  linear  program  that  was  proposed  in  [12]. 

This  method  is  based  on  the  fact  [14]  that  the  linear  program  (17)  is 

solvable  if  and  only  if  the  quadratic  program 

minimize  f-x^x  +  cTx  subject  to  Ax  >  b,  x  >  0  (28) 

xcR"  2 

is  solvable  for  all  ee  (0,e)  for  some  e  >  0.  Furthermore  the  unique 
solution  of  (28)  is  independent  of  e  for  ee  (0,e)  and  is  the 

closest  solution  of  the  linear  program  (17)  to  the  origin  in  the  2-norm 

[14].  Note  that  e  may  be  infinite  in  some  special  cases.  Problem 
(28)  can  be  solved  then  by  the  QPSOR  algorithm  of  Section  3.  From  (8) 
the  dual  to  the  quadratic  program  (28)  is 

maximize  -l-x^x  +  b^u  subject  to  ex  -  A^u  -  v  +  c  *  0,  (u,v)>0  (29) 

(x,u,v)€Rn+m+n  c 

which  upon  elimination  of  x  by  using  the  constraint  relation 

x  *  ~(ATu+v-c)  (30) 


-15- 


gives 

IT  2  T 

minimize  4- 1| A  u+v-c II  -  eb  u  subject  to  (u,v)  >  0 
(u,v)cR”+"  2  ■  "  " 

Note  that  (31)  is  the  classical  exterior  penalty  function  [7] 
associated  with  the  dual  linear  program  (18).  However  the  perturba¬ 
tion  results  of  [14]  give  the  stronger  result  that  c  in  (31)  need 
not  approach  zero  in  order  for  x  defined  by  (30)  to  be  a  solution 
of  (17).  In  other  words  if  we  let  (u(e),  v(e))  be  a  solution  of 
(31)  for  £  €.  (0,e)  then  x  *  ^(ATu(e)  +  v(e)  -  c)  is  independent  of 
e  and  is  the  closest  solution  of  the  linear  program  (17)  to  the 
origin  in  the  2-norm.  Note  however  (u(e),  v(e))  need  not  be  a 
solution  of  the  dual  linear  program  (18)  for  ee  (0,c)f  but  each 
accumulation  point  of  {(u(e- ),  v(e. )))  will  be  a  solution  of  (18) 
if  {Cj>  is  a  decreasing  sequence  converging  to  zero.  We  can  now 
solve  (31)  by  a  sparsity-preserving  algorithm  which  follows  directly 
from  the  QPSOR  algorithm  of  Section  3  by  replacing  D  by  el. 

LPS0R2  Algorithm 

Choose  (u°,v°) €  Rj*",  we  (0,2)  and  e  >  0.  Having  (u*,v*) 
determine  u^+1,  vi+^  as  follows: 


X 

for  j>l 

vi+1  a  (v^aV+V-O) 


U*+W - “ 2(V  l  ?  (AJ).  uW-c)-eb 

j  J  IJA.II2  j  1=1  1  1  £*j  11  J 

J  for  .  . 


)) 


jcl  9  •  •  •  •  |H) 


(3)) 


(32) 
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Note  that  this  LPS0R2  algorithm,  unlike  the  algorithms  proposed 
in  [12],  will  preserve  any  sparsity  the  matrix  A  may  have  and  there 
is  no  need  to  compute  AA^  as  was  done  in  [12]  and  thereby  destroying 
any  sparsity  that  A  may  have  had. 

The  following  convergence  theorem  follows  directly  from  the  con- 
*  vergence  theorem  of  the  QPSOR  algorithm,  Theorem  2  and  the  perturbation 
results  of  [14]. 

Theorem  4:  LPS0R2  Convergence 

(i)  Let  the  linear  program  (17)  have  a  solution.  There  exists  a 

real  positive  number  e  such  that  for  each  e  in  the 

interval  (0,e),  each  accumulation  point  (u,v)  of  the 
•  • 

sequence  { (u1 , v1 ) >  generated  by  the  LPS0R2  algorithm  (32) 
solves  (31)  and  the  corresponding  x  determined  by  (30)  is 
Independent  of  e  and  is  the  (unique)  solution  of  the  linear 
program  (17)  which  is  closest  to  the  origin  in  the  2-norm, 
(ii)  If  in  addition  to  the  assumptions  of  part  (i)  the  constraints 
of  the  linear  program  (17)  satisfy  the  Slater  constraint 
qualification  (15)  then  the  sequence  {(u^v1)}  of  the 
LPS0R2  algorithm  (32)  is  bounded  and  has  an  accumulation 
point  for  each  ee(0,e). 
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